    %% Calculate the Implied Tax on Debt for states in which the constraint does not bind. Set the rest to 0
Tau1p=ones(NBG,NB,NSS,NBG); %THis is Tau + 1
%Tau1p_alt=ones(NBG,NB,NSS,NBG); %THis is Tau + 1
price=price_cc;
bmax=2;
bmin=-max(BGgrid);

% POLICIES IN REPAYMENT
totalc=(omega*c_big.^(-ita)+(1-omega)*yn_big.^(-ita)); 
mupx=real(omega*totalc.^(sigma/ita-1/ita-1).*(c_big.^(-ita-1)));

                        bpmax=(1./q_big).*KAPPAS_big.*(price.*yn_big + yt_big);
                        bpmax(bpmax>bmax)=bmax;
                        bpmax(bpmax<bmin)=bmin;
                        cbind=yt_big+q_big.*bpmax-(1-PVD_big).*b_big+Transfer;
%Policies in Default
totalc_bad=(omega*c_def_2d.^(-ita)+(1-omega)*yN2D.^(-ita)); 
mup_badstanding=(omega*totalc_bad.^(sigma/ita-1/ita-1).*(c_def_2d.^(-ita-1)));

                        price_bad=((1-omega)/omega*(c_def_2d./yN2D).^(1+ita));
                        bpmax_bad=(1./q2D).*KAPPAS2D.*(price_bad.*yN2D + yT2D);
                        bpmax_bad(bpmax_bad>bmax)=bmax;
                        bpmax_bad(bpmax_bad<bmin)=bmin;

                       
    %% COMPUTE EMU
      Premu=Def_func.*reshape(repmat(mup_badstanding(:),1,NBG)',NBG,NB,NSS)+Repay.*sum(mupx.*G_big,4);                                                                            
      Premu_bad=(1-theta)*mup_badstanding+theta*squeeze(Premu(zeroBG_ind,:,:)); % YOU REE ENTER WITH NO DEBT
               

             parfor j=1:NSS
              Ptemp=Prob(j,:)';
              bc_bad=squeeze(bp_def_2d(:,j)); 
              bc_tot_bad=bc_bad(:);
              emu_temp_bad=PVDm1_2d.*Premu_bad*Ptemp;
              emu_interp_bad=griddedInterpolant(bgrid1d,emu_temp_bad,'linear','linear');
              emu_mat_bad(:,j)=emu_interp_bad(bc_tot_bad);
                          for inext=1:NBG                                      
                          bc=squeeze(bp(:,:,j,inext)); 
                          bc_tot=bc(:);
                          emu_temp=squeeze(PVDm1_2d.*squeeze(Premu(inext,:,:)))*Ptemp;
                          emu_interp=griddedInterpolant(bgrid1d,emu_temp,'linear','linear');  
                          emu_matR(:,j,inext)=(emu_interp(bc_tot));
                          end
                end
                
           emu_mat=reshape(emu_matR(:),NBG,NB,NSS,NBG);   


%% Calculate Tax settting Tax rate at zero when the constraint binds
Tau1p=(q_big.*mupx)./(beta*emu_mat);

                parfor BGi=1:NBG %BGi here is the index of the government debt
                        for i=1:NB 
                            for j=1:NSS
                                for inext=1:NBG
                      
                                       if bpmax(BGi,i,j,inext)-bp(BGi,i,j,inext)<1e-4
                                        Tau1p(BGi,i,j,inext)=1;
                                     end
                                end
                            end
                        end
                end


Tau1p(Tau1p<1)=1; 

 Tau1p=Tau1p-1.;


%% THERE IS ALSO THE TAU in BAD STANDIN
Tau_bad=q2D.*mup_badstanding./emu_mat_bad;
              for i=1:NB 
                            for j=1:NSS
                                
                                       if bpmax_bad(i,j)-bp_def_2d(i,j)<1e-4
                                        Tau_bad(i,j)=1;
                                     end
                                end
                    end
 max(Tau_bad(:))
min(Tau_bad(:))
Tau_bad(Tau_bad<1)=1; 
Tau_bad=Tau_bad-1.;
